function kl = gausskl(m1, s1, m2, s2)

kl = log(det(s2)/det(s1)) + trace(inv(s2)*s1) + (m1-m2)*inv(s2)*(m1-m2)' -length(m1);
kl = kl * 0.5;
